Exponential decay properties of Wannier functions and related quantities 
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The spatial decay properties of Wannier functions and related quantities have been investigated 
using analytical and numerical methods. We find that the form of the decay is a power law times 
an exponential, with a particular power-law exponent that is universal for each kind of quantity. 
In one dimension we find an exponent of —3/4 for Wannier functions, —1/2 for the density matrix 
and for energy matrix elements, and —1/2 or —3/2 for different constructions of non-orthonormal 
Wannier-like functions. 

PACS: 71.15.Ap, 71.20.-b, 71.15.-m 
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A growing interest in localized real-space descriptions 
of the electronic structure of solids has been motivated 
by the development of computationally efficient "linear- 
scaling" algorithms and by the desirability of a local 
real-space mapping of chemical [^,^ and dielectric [^] 
properties. A primary avenue to such a description is 
the use of Wannier functions |@-^ (WFs), i.e., a set of 
localized wavefunctions WR(r) obtained from the Bloch 
functions ipkir) by a Fourier- like unitary transformation. 
A closely related approach is to represent the electronic 
structure in terms of the density matrix n(r, r'). It is 
thus not surprising to find considerable recent interest 
in the localization properties of the WFs Q and of the 
density matrix |lO|JTl| ]. 

In a classic 1959 paper, Kohn proved, for the case of a 
centrosymmetric crystal in one dimension (ID), that the 
WFs have an "exponential decay" w{x) « e"''^, where h 
is the distance of a branch point from the real axis in the 
complex- A: plane H. More precisely. 



lim w{x) e'^ 



0, 



q < h 
q> h 



(1) 



The density matrix has a similar decay n{x,x') 



The exponential decay of the WFs has since 



-h\x—x'\ 

been proven for the general ID |1^ and single-band 3D 
p3[ cases, and that of the density matrix (more pre- 
cisely, of the band projection operator) has been proven 
in general jl^. The energy matrix elements E{R) — 
{wr\H\wo), with WfiXx) = w{x — R) and R = la a lat- 
tice vector, are also expected to have a similar decay, 
E{R) - e~^". 

The purpose of this Letter is to address two questions. 
First, Eq. (Q) allows considerable freedom; in fact, it is 
consistent with 



w{x) 
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(2) 



for any exponent a, i.e., a decay which could be faster 
{a > 0) or slower [a < 0) than pure exponential. Does 
such a power-law prefactor exist, and if so, what is the 
exponent a? Second, it has long been understood that 
relaxation of the orthogonality constraint {'Wq\'Wii) ^ Sq u 
can give "more localized" Wannier-like functions [^-p^ . 



In what sense are these more localized - a larger h, or 
a larger a for the same h, or only a smaller prefactor 
of the tail? We show that the power-law prefactors of 
Eq. (H) do exist, and that the various quantities have 
a common inverse decay length h but different expo- 
nents a. In ID we find that a — 3/4 for usual (or- 
thonormal) WFs, a = 1/2 for n{x,x') and E{R), and 
a = 1/2 or a = 3/2 for two different constructions of non- 
orthonormal Wannier-like functions (NWFs) . The NWFs 
of superior decay ('^ x~^/^e~''^) can be constructed by 
a projection method as duals to a set of trial functions. 
These results may have important implications for the de- 
sign and implementation of efficient real-space electronic- 
structure algorithms. 

We first review the central results of the pioneering 
work of Kohn Q , who considers a centrosymmetric po- 
tential of period a in ID. The WFs are constructed as 



Wn{x - R) ^ WnB.{x) = 
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ipnk{x)dk (3) 



with the phases of the Bloch functions ipnk chosen as in 
Sec. 6 of Rcf . [|j . The exponential decay of the WFs is 
then governed by the positions of branch points in the 
"complex band structure" En{k) constructed by regard- 
ing complex En to be a function of coniplex k via ana- 
lytic continuation from the real axis |^,^. Specifically, 
there is a Riemann sheet En{k) for each band n, and the 
branch points fc„ are the points at which the sheets are 
connected, i?„(fc„) = i?„+i(fc„). These are located at 



7r/a ± ih 
±ih. 



„, n even 
n odd 
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and at translational image locations kn"' = fc„ -t- 2'Km/a 
for integer m. En[k) and V'ra(fc) are thus analytic func- 
tions in the strip |Im(fc)| < hn where /i„=min(/i„_i, /i„). 
Kohn's main result |^ is that the decay of the WF for the 
n'th band is as Wn{x) ~ e~''"l^l in the sense of Eq. (|l). In 
what follows we restrict our attention to the bottom band 
(n=0), for which = ho (henceforth just h). The rele- 
vant branch point in the upper half-plane is kg = TT/a-\-ih 
and the expected Wannier decay is w{x) w 
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FIG. 1. Decay of normalized WFs w{x) and NWFs v{x) 
and y{x). (a) ID model (see text) with b=0.3 and Vo=— 10. 
(b) 3D Si plotted along the [110] direction. 



To confirm this decay numerically, we first choose a 
simple ID model Hamiltonian having a periodic poten- 
tial U{x) — X]m^at(a; — TTT-a) constructcd as a sum of 
Gaussian "atomic" potentials Vi,t{x) — {Vo/b^/Tr)e~^^/''^ . 
Here a is the lattice constant and Vq and b control the 
depth and width of T4t. We choose units such that 
m=h=e=l and keep a=l and Vo=~10 fixed while ad- 
justing 6 to vary the gap. The Bloch functions are com- 
puted on a mesh of 200 k points by expanding in 401 
plane waves and the WF at i? = is then constructed 
according to Eq. (^) using 128-bit arithmetic. 

The resulting decayof the WF for 5=0.3 is shown as 
the solid line in Fig. |l|(a). In this semilog plot, the ap- 
proximate linearity of the peaks is consistent with the 
expected exponential decay, but there is a slight curva- 
ture that can be analyzed further. To do so, we first 
computed the En{k) along it /a + in for real k and de- 
fined h to be the value of k at which Eq = Ei. For 6=0.3 
we find ft,= l. 28869. In Fig. ||(a) we then plot (diamonds) 
hx + ln \w(x)\ vs. ln(a;) for each peak of In |w(x)|. A pure 
exponential decay w{x) « e~^^ should yield a horizontal 
line in such a plot; instead, the data appears linear with 
a slope of —3/4, indicating that 



w{x) 



A similar plot (not shown) for 



E{R) = {wr\H\wo) 



a 

2^ 



dk e'^''E{k) 



(5) 



(6) 



suggests that E{R) shares the same inverse decay length 
h but has a different power-law exponent. 
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FIG. 2. (a) As in Fie. |l|(a) but plotted so that slope reveals 
exponent —a of Eq. (^j. (b) Same for fe=0.6 (nearly-free elec- 
tron case) showing crossover. Pluses, diamonds, and squares 
represent v, w, and y, respectively. 



Naturally h changes if the potential parameter b is varied, 
but we find that the power-law exponents of —3/4 and 
—3/2 do not. It thus appears that these exponents are a 
universal feature of electron bandstructures in ID. 

In order to gain an analytic understanding of this be- 
havior, we consider first the simpler case of the energy- 
band Fourier transform E{R) ^ E{k). Kohn showed 
that the expansion of E{k) about ko = ir/a + ih takes 
the form g 



Eik)=Eo+-/ik-ka)'/^ + 



(8) 



with higher terms of order (fc — fco)^, {k — fco)"^/^, etc. 
The form of this expansion arises from the requirement 
that E{k) come back to itself if k traverses a closed path 
winding twice around /cq, consistent with the picture of 
two Riemann sheets touching at kg. 

Now there are well-known mathematical results that 
relate the behavior of a function near a branch point to 
the asymptotic decay of its Fourier transform ||l^. The 
following lemma is useful here. Let /(fc) be a periodic 
function f{k) = f{k + 27r/a) that has a leading behavior 



I{k) = /o + 7 V{k - fco)]'' 

when expanded at the branch point fco = tt/o + ih. 
Fourier series coefficients are given by 



(9) 



Its 



F{x) = / f{k) 6*'=^ dk 
J Co 



(10) 



E{R) w 



-hR 
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at a; = ma for integer m. As shown in Fig. ^, the contour 
Co initially lies along the real axis. However, f{k) e'*^^ 
is invariant under k ^ k -\- 27r/a, and assuming that no 
other branch points or poles intervene, the contour can 
be deformed to become Ci as shown in Fig. ||. The ex- 
ponential smallness of e*'^'^ for large x kills the integrals 
along the horizontal segments, and the dominant con- 
tribution to the Ci integral comes from the vicinity of 
/cq. Using the contour-integral definition of the Gamma 
function [ITsI, 
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FIG. 3. Branch points (x), cuts (dashed Unes), and inte- 
gration contours (Co and Ci) in the complex-fc plane. 

|F(a;)| ~7B(/3)x-(i+«e-''" , (11) 

where Bg = 2 sin(/37r) r(l + 

Eq. (^l]) now allows us to understand the observed be- 
havior of quantities such as E{R) and w{x). For ex- 
ample, since E{k) in Eq. (||) has /3 = 1/2, we confirm 
that E{R) « R-"e-''" with a = 1 + /3 = 3/2. Simi- 
larly, to understand the decay of w(x), we need to know 
the behavior of tpk{x) regarded as a function of k near 
the branch point feg. Once again Kohn Q provides the 
needed result V'fc ~ (fc — ^o)~^^'*- Sure enough, (3 — —1/4, 
gives a decay w{x -I- i?) ~ R~^/^e~'^^ for small x and 
R ^ a. In other words, ?«(x) « x~'^/^e~''^ for large x, 
as obtained numerically from Fig. ^(a). 

We can summarize the information about both the k- 
and x-dependence of 4'k{x) near the branch point as 

M^) = M^) 9"'/' + ^i(^) q^'^ + ■■■. (12) 

where q — i{k — ko). (All such terms have powers 
that are odd-integer multiples of 1/4, consistent with 
ipkix) —ipkix) when traversing a closed path wind- 
ing twice around fcg.) Aq{x) and Ai{x) are real functions 
obeying An{x + a) = e*'""" An{x). 

The locality of the density matrix (i.e., the band pro- 
jection operator) is also very important. For example, 
many linear-scaling algorithms are based on a direct so- 
lution for the density matrix [p]j2|p9|] . We can write 

"-(2;',2;) = TT" / i^-k{x')i^k{x)dk (13) 
2^ J Co 

where, following Kohn we have substituted ij^t [x') 
by Tp^kix') in order that the integrand of Eq. ( [l^ ) 
should remain analytic off the real axis. The behavior of 
ijj^k{x) near the branch point is 'ijj-^kix) ~ Aq{—x) [i{k — 
ko)]^^''^. The integrand of Eq. ( [13|) then takes the form 
i:-k(y)Mx) « Ao{-x')Ao{x)[iik - ko)]-^^^. Applying 
Eq. ( [ll| ) yields n(0, x) « x~^^'^ e~^^ for large x, and more 
generally, n{x',x) ^ (x — x')"^/^ e~''(^~^ ^ for a; ^ x' . 
This a = 1/2 behavior of the decay has been confirmed 
from numerical plots (not shown) similar to Fig. ^a). 

We have so far shown that E{x), w{x), and n{0,x) all 
have a decay of the form x~°' e~^'^ with a common h but 
with different (universal) exponents a£;=3/2, a^=3/4, 
and q;„=1/2. The energy matrix elements thus have the 
fastest decay, and the density matrix the slowest. 



One may next ask whether it is possible to find 
non-orthonormal Wannier-like functions (NWFs) with a 
faster decay than those of the orthonormal WFs w{x) 
[p^p^ . We explore this question in the context of band- 
projection methods po|-p^. We find that a naive ap- 
plication of the projection technique actually generates 
NWFs with a slower decay, while a modified "dual con- 
struction" approach does give improvement as measured 
by the exponent a. 

The basic idea of the projection technique is to start 
with a trial function t{x) and generate a Wannier-like 
function v(x) by acting with the band-projection op- 
erator P = Efc IV'fc)(V'/c|, i-e., \vR)=P\tR). Here \tR) 
corresponds to the translational image t{x — R) of t{x) 
in cell R = na, and similarly for |u_r,). The trial func- 
tions can be Gaussian functions, atomic or molecular or- 
bitals, etc. The \vi{) are NWFs having overlap Sqb. = 
{v{)\vr) = (io| P\tR). Numerical investigations on C and 
Si by Stephan and Drabold indicated that the projected 
functions v{x) are not more localized than the true WFs 
'w{x) [Q. This should not be surprising; introduction of 
NWFs may give flexibility to generate more localized or- 
bitals, but this flexibility needs to be used to advantage. 
To do so, we introduce dual functions y{x) defined via 
\vo)=Y.r.{S'^)or\vr), so that (yoki?,) = 5f)R and also 
(2/0 K-r) = i^O-R- This latter equation means that y{x) is 
orthogonal to the trial function at every site except i?=0, 
suggesting that y{x) may be especially well localized. 

Numerical tests of the decay of (normalized versions 
of) v{x) and y{x) are shown as dashed and dotted curves 
respectively in Figs, ^(a) and |2|(a). The trail function 
used is a (5-function on the atomic site, but use of other 
narrow trial functions gives similar results. It clearly 
appears that a = 1/2 and 3/2 for v{x) and y{x) respec- 
tively, to be compared with a = 3/4 for 'w{x). Thus, the 
simple projected functions v{x) actually have a slower de- 
cay that the WFs ■u'(a;), but the duals y{x) have a much 
faster decay than either of them. 

These results can be explained by the complex analy- 
sis of the Bloch-like functions Vk{x) and yk{x) that are 
related to v{x) and y{x) in the same way that tpkix) is 
related to w{x). Defining 

/oc 
ij-k{x)t{x)dx , (14) 
-00 

it follows from \vk) = \ipk){'ipk\t) that \vk) = rik\->pk)- 
Also the Fourier transform of S{0,R) can be seen to be 
S{k) = rj'^ik), so that \yk) — {ipk) /'ri{k). In the vicinity 
of fco we have 

f]{k) = 770 [i{k ~ ko)]-^^^ + m Hk - fco)]'/' + • • • (15) 

where rjn = An{—x) t{x) dx. Moreover, 

Vk{x) = r]oAo{x)q^'^/^ H , 

Wk{x) ^ Ao{x) q-^/^ + ■ ■ ■ , 

yk{x) = - {Aoix)+Mx)q'/' + ---] , (16) 
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where q = i{k — ko) and Ai{x) = Ai{x) — (?7i/??o) Ao{x). 
The leading term in yk{x) gives no singularity, so the 
real-space decay is determined by the behavior of the 
next term for which a = P + 1 = 3/2. To be explicit, we 
can define An{x) = An{x) e^'^ so that A is anti-periodic, 
A{x + a) = —A{x), and for large x we find 



w{x) ~ B_i/4Ao{x)x~^/^ : 
v{x) ~ B_i/2'noAa{x)x 



-hx 



-1/2, 



y{x) ~ {Bi,2h^)Ai{x)x ^'"^e 



hx 



(17) 



The above conclusions regarding y{x) rely on the ab- 
sence of zeros of r]{k) inside the strip —h < Im(fc) < h. If 
such zeros exist, ykix) = ipkix)/'r]k may have new singu- 
larities and y{x) will then have poor decay compared to 
other NWFs. We find that this problem does not arise 
when using t{x) = d{x) or a narrow Gaussian, but can 
be triggered by use of a too-wide Gaussian for t. 

In view of n(x',x) = J^i'^ii^') '^ii^) it may appear 
surprising that n{x',x) decays more slowly than wlx) 
{an — 1/2 vs. au,=3/4). The representation of n via w 
thus has some advantages. Better yet, perhaps, one can 
represent n(x' ,x) — ^ij Vii^') Vji^)- Here the slow 
decay has been transferred to a simple matrix quantity 
(as=l/2) but the NWFs decay very quickly (aj,=3/2). 

Is it possible to find a NWF with an even faster decay 
than a;-^/^e~'"? Yes; define a new NWF Zk ~ f{k)yk 
where f{k) is analytic in the strip |Im(fc)| < h and has 
simple zeros at the branch points (2n + 1)tt ± ih. The 
function f{k) = l-l-cos(fca)/ cosh(ft,a) is a good candidate 
p3| . Then the leading singularity of Zk is as {k — fco)^^^, 



and we expect z{x) 



5/2g hx ^ have confirmed 



numerically that this works. However, since the multi- 
plication by f{k) in fc-space corresponds to a convolution 
in real space, the resulting z{x) is actually broader than 
y{x) or w{x) by almost any other measure (e.g., second 
moments [^). Thus, this strategy may be counterpro- 
ductive in practice. 

Before leaving the ID case, we make two brief com- 
ments. First, the extension to the case of non- 
centrosymmetric potentials in ID is not difficult, and 
the results (including values of the a exponents) are un- 
changed. Second, there is an apparent paradox concern- 
ing the nearly- free electron limit. For free electrons the 
occupied portion of the band gives w{x) ^ sm.{kFx)/kFX, 



One may expect this to go over to ' 



-hx 



in the nearly-free case, but this would be inconsistent 
with our general result a = 3/4. Actually we find there 
is a crossover behavior, as shown in Fig. 0(b), with a—1 
behavior for x << Xc and a=3/4 in the true large- a; tail. 
The crossover distance Xc increases as the gap decreases; 
as the gap closes, Xc ^ oo and h ^ 0. 

Before concluding, we briefly discuss the 3D case. Here 
the 6-dimensional space of complex {k^, ky,kz) makes the 
formal analysis difficult. We have carried out a numerical 
calculation of WFs and NWFs for Si using an empirical- 
pseudopotential scheme starting from four bond-centered 
trial functions. The results are plotted in Fig. 0(b). Plots 



of /ix -|- In |_F(x)| vs. ln(a;) (not shown) again show linear 
behavior, with slopes that appear consistent with the ID 
values of a = 3/4, 1/2, and 3/2 for w, v, and y, respec- 
tively. However, in this case we cannot afford to go to 
very large x values, and we suspect that there may be a 
crossover to larger a values in the far tails. We leave this 
as a question for future investigations. 

To conclude, we find that in ID the asymptotic behav- 
ior of WFs and related quantities can all be expressed 
as x~'^e~'^^ with a common h, and with exponents a 
that take on universal rational values depending on the 
type of singularity of the relevant function at the branch 
points in the complex-fc plane. It is surprising that this 
behavior has gone unnoticed since Kohn's seminal 1959 
paper. The consequences for linear-scaling calculations, 
and localized real-space representations of electron struc- 
ture more generally, remain to be fully explored. 
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